Direct measurement of the 3He+ magnetic moments

Helium-3 has nowadays become one of the most important candidates for studies in fundamental physics1–3, nuclear and atomic structure4,5, magnetometry and metrology6, as well as chemistry and medicine7,8. In particular, 3He nuclear magnetic resonance (NMR) probes have been proposed as a new standard for absolute magnetometry6,9. This requires a high-accuracy value for the 3He nuclear magnetic moment, which, however, has so far been determined only indirectly and with a relative precision of 12 parts per billon10,11. Here we investigate the 3He+ ground-state hyperfine structure in a Penning trap to directly measure the nuclear g-factor of 3He+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{I}^{{\prime} }=-\,4.2550996069(30{)}_{{\rm{stat}}}(17{)}_{{\rm{sys}}}$$\end{document}gI′=−4.2550996069(30)stat(17)sys, the zero-field hyperfine splitting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{{\rm{HFS}}}^{\exp }=-\,8,\,665,\,649,\,865.77{(26)}_{{\rm{stat}}}{(1)}_{{\rm{sys}}}$$\end{document}EHFSexp=−8,665,649,865.77(26)stat(1)sys Hz and the bound electron g-factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{e}^{{\rm{\exp }}}=-\,2.00217741579(34{)}_{{\rm{stat}}}(30{)}_{{\rm{sys}}}$$\end{document}geexp=−2.00217741579(34)stat(30)sys. The latter is consistent with our theoretical value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{e}^{{\rm{theo}}}=-\,2.00217741625223(39)$$\end{document}getheo=−2.00217741625223(39) based on parameters and fundamental constants from ref. 12. Our measured value for the 3He+ nuclear g-factor enables determination of the g-factor of the bare nucleus \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{I}=-\,4.2552506997(30{)}_{{\rm{stat}}}(17{)}_{{\rm{sys}}}(1{)}_{{\rm{theo}}}$$\end{document}gI=−4.2552506997(30)stat(17)sys(1)theo via our accurate calculation of the diamagnetic shielding constant13 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sigma }_{{}^{3}{\mathrm{He}}^{+}}=0.00003550738(3)$$\end{document}σ3He+=0.00003550738(3). This constitutes a direct calibration for 3He NMR probes and an improvement of the precision by one order of magnitude compared to previous indirect results. The measured zero-field hyperfine splitting improves the precision by two orders of magnitude compared to the previous most precise value14 and enables us to determine the Zemach radius15 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{Z}=2.608(24)$$\end{document}rZ=2.608(24) fm.


Supplementary Information
1 Supplement -experiment

Magnetic inhomogeneity
The quadratic inhomogeneity of the magnetic field in the PT B 2,P T was determined by exciting the modified cyclotron mode with a resonant signal applied to a split trap electrode and measuring the shift of the axial frequency as well as the shift of the modified cyclotron frequency The latter includes besides the term proportional to B 2,P T an additional relativistic contribution, so that B 2,P T = 1.17(4) T m −2 can be extracted from the ratio of both shifts (Fig. S1).

Electrostatic anharmonicity
Optimizing the tuning ratio TR = Vc Vr , defined as the ratio of the voltages applied to the center electrode V r and the neighbouring correction electrodes V c , allows to eliminate anharmonic contributions C 4 and C 6 to the electrostatic potential Φ(z) = k C 2k z 2k , which would cause energy dependent shifts of the eigenfrequencies [1]. For this purpose the magnetron mode is excited to different energies E − for different tuning ratios and the resulting shift of the axial is measured. By this method the anharmonic contributions were minimized to C 4 = 0(5) · 10 4 m −4 and C 6 = 2(1) · 10 −6 mm −6 . The axial frequency shift in Eq. (S3) includes besides the electrostatic anharmonicities also a B 2 dependent term, which, however, is negligible assuming the measured value of B 2,P T . As a cross-check, the tuning ratio was additionally optimized by heating the axial mode using white noise and measuring the resulting shift of the axial frequency, which is B 2 independent δν z ν z = 3 4 The optimal tuning ratios found by both methods agree within 1σ.

Axial temperature
In order to determine the axial temperature T z in the PT, the axial and cyclotron modes were coupled in the PT at the sideband. As a result, the cyclotron mode thermalizes with the axial detection system at T + = T z ν+ νz . The ion is afterwards transported to the AT where the axial frequency is measured. Due to the magnetic bottle in the AT, the axial frequency is shifted proportional to the cyclotron energy E + = k B T + , see Eq. (S1). A maximum likelihood fit of a Boltzmann-distribution to the axial frequencies therefore gives the temperature T + and thereby T z = 11.7(7) K. As a cross-check, the axial temperature in the PT was additionally determined by exciting the modified cyclotron mode and afterwards measuring the axial frequency in the PT. After repeating this 200 times, the temperature can be calculated from the standard deviation and average of the axial frequency distribution as described in [2]. This method requires inserting D 4 = dC 4 /dTR as obtained theoretically from the trap geometry D 4 = −3.5 · 10 9 m −4 and the measured B 2,P T and gives a consistent result for the axial temperature T z = 12(1) K,

Time standard
All frequency generators and FFT signal analyzers are locked with a Stanford Research Systems FS725 Rubidium 10 MHz frequency standard. As the determination of E HFS requires an absolute frequency measurement, the Rubidium standard is itself locked to a GPS pps signal.

Lineshape
Due to the oscillation of the ion in the residual magnetic bottle in the PT B 2,P T , the resonance frequency varies during the spin excitation time T exc = 120 s. Accordingly, the linewidth parameter δω(B 2,P T , T z ) characterizing the spread of the resonance frequency is dependent on the magnetic bottle and the average axial amplitude ⟨z 2 ⟩ ∝ T z .
Here, Ω is the Rabi frequency and γ is the coupling of the ion to the axial resonator. The correspondent spin-flip probability P 1/2 (∆) as a function of the detuning ∆ = ν RF − ν L of the applied drive frequency ν RF from the resonance frequency ν L is then The expected resonance frequency ν L (B) follows from the double-dip measurement of the cyclotron frequency, which is normally distributed with a relative width of 1.7 · 10 −9 . The underlying function describing the spin-flip probability is therefore convoluted with the Gaussian of the magnetic field measurement. For small Rabi frequency Ω or small B 2,P T T z , the jitter of the magnetic field measurement dominates and the spin-flip probability can be approximated as a Gaussian. In order to estimate the systematic uncertainty due to the unknown deviation from a Gaussian profile, the two possible spin-flip probabilities P 1/2 (∆) convoluted with a Gaussian are fitted with a Gaussian line to extract the shift of the center-frequency and the correspondent difference of the deduced g factors and E HFS (Fig. S5).
The fitted Rabi frequencies are different for both lineshapes and shown for the electronic transitions |1⟩ ↔ |3⟩ and |2⟩ ↔ |4⟩ in Fig. S6 and for the nuclear transitions |1⟩ ↔ |2⟩ and |3⟩ ↔ |4⟩ in Tab. S1. Each measurement is ascribed a systematic uncertainty given by the maximum of the differences between a Gaussian fit and a fit with Thus, the results for the electronic and nuclear parameters include 1540 and 490 data points, respectively, with two data points measured per hour in both cases.
The coupling of the axial and modified cyclotron modes during the double-dip measurement leads to a modulation of the axial oscillation amplitude and thus of δω with the Rabi frequency of the coupling Ω DD = 22 Hz, causing sidebands in the resonance lineshape. A significant shift due to this effect was excluded by measuring the |2⟩ ↔ |4⟩ resonance once with and without a double-dip during the spin-flip drive at the same microwave power, giving results that agree within 1σ. For this purpose the cyclotron frequency measurement during the spin-flip drive was replaced by the average of the measured frequencies before and after applying the drive, as these values were consistent for all measurement runs. a b Figure S4: Results for the shielded nuclear g factor and hyperfine splitting with a maximum likelihood fit of a Gaussian, combining one resonance of each nuclear transition with microwave powers P 12 and P 34 . As the electron g-factor is known precisely by theory, for the final result g e is fixed to the theoretical value in the fits of the nuclear resonances. Plotted are the results if g e = g e,theo is assumed and when g e is changed by 3σ of the theoretical value, i.e. g e /g e,theo − 1 = −6 · 10 −13 for all microwave powers. In order to make the overlapping errorbars better visible, the datapoints are depicted slightly shifted along the y-axis. Changing g e by 3σ shifts the fitted g ′ I and E HFS by two orders of magnitude less than their statistical uncertainties. Table S1: Fit results for the Rabi frequencies Ω 1 and Ω 2 of the two measured nuclear transitions assuming lineshapes χ 1 and χ 2 , respectively, which correspond to a systematic uncertainty shown in Fig. S5. As g I and E HFS are determined from a combined fit of one resonance of each nuclear transition, the combined uncertainty of both transitions are given for these parameters. The g e factor for 3 He + was calculated similarly to that of 4 He + in Ref. [5]. The contributing terms are shown in Table S2. The most notable difference between the two isotopes is the recoil contribution due to a different nuclear mass. Calculating the nuclear mass using the GRASP code [6], with the He-3 mass from the recent atomic mass evaluation [7] as input parameter, resulted in recoil corrections which are identical to all digits specified to results obtained with the helion mass from [8]. The result for the leading recoil correction to O (Zα) 2 [9][10][11] differs from the all-order in Zα result [12] given in Table S2 by 1.18 × 10 −12 . For the calculation of the higher-order in m M recoil corrections, we used a formula which is all-order in the mass ratio [10] (later confirmed in [13,14]). We note c d e f Figure S5: Difference of the fitted electron g-factor, nuclear g-factor and E HFS with the models P 1 (A,C,E) and P 2 (B,D,F) compared to a Gaussian as a function of the Rabi frequency extracted from the fit. Gaussian P2 g e /g e,theo -1 here that the result up to O m M 2 [9,11] deviates by only 1 × 10 −14 from the all-order result. Apart from that, several recent results were taken into account. Free-QED contributions have been reevaluated, namely, the fourloop coefficient was evaluated to high precision in Ref. [15], and the five-loop term has been updated in Ref. [16].
Recently computed two-loop binding contributions of O (Zα) 5 [17,18] are very small for Z = 2. Higher-order corrections to two-loop magnetic loop vacuum polarization diagrams [19] turned out to be negligibly small for Z = 2.
Binding corrections of order (Zα) 2 were calculated using a simple formula from Ref. [20], which is a confirmation of earlier derivations of this quantity [9][10][11]. The uncertainty due to uncalculated binding corrections of order (Zα) 4 and higher to Feynman diagrams with at least three loops was estimated by adapting a method from Ref. [21].
The uncertainty due to uncalculated binding corrections of order (Zα) 6 and higher to Feynman diagrams with at two loops was calculated using the methods from [21] and [17], with the larger of the two estimates given as the uncertainty of this contribution in Table 2. The FS correction was calculated nuclear model independently using analytic formulas from [22][23][24] and the rms charge radius from [25]. The muonic vacuum polarization correction [26] as well as our estimates for the nuclear deformation [27], nuclear susceptibility [28] and our calculated nuclear polarization value are also much smaller than 1 in the last digit shown in Table S2. Finally, the latest all-order evaluation of the one-loop self-energy contribution has been included [29]. With this, the uncertainty of the total g e factor is dominated by the uncertainty of α.
It is worth pointing out that an independent calculation of the 5-loop contribution to the free electron anomaly [30] is in disagreement with the result from [8] by four times its uncertainty as specified in [8].
Our calculation of the bound-electron g factor relies on the latest CODATA value for the fine-structure constant.
It is worth mentioning that since the CODATA 2018 compilation, α has been measured with a better accuracy than specified in CODATA 2018 [31]. This new value of α disagrees with the CODATA value by 5 σ. However, since the uncertainty due to the (CODATA) fine structure constant is three orders of magnitude smaller than the experimental uncertainty of the work presented here, we do not expect this deviation to be relevant for this work.

Theory of the hyperfine splitting
The energy difference between the hyperfine levels of the ground state, with the total angular momenta F = 0 and F = 1 (electronic and nuclear spins oriented anti-parallel and parallel, respectively), is  with (see Ref. [37,49]) with the quantities defined in the manuscript. This leads to with Here, I = 1 2 is the nuclear spin [50], and the angular momentum of the electron is j = 1 2 . The relativistic factor for the 1s state is with γ = 1 − (Zα) 2 , and the mass factor is M = 1 + me , with M N being the nuclear mass.
An approach to express the finite-size correction δ FS to the HFS for low Z was developed in Refs. [51,52], with δ E and δ M corresponding to the finite size charge and magnetization distribution corrections, respectively. The formula for the combined effect, as derived in Refs. [51,52], is as follows: The sum of the leading nonrelativistic corrections −δ nr E − δ nr M is then identified as the Zemach correction [53] The Zemach radius is defined as with the electric G E (q 2 ) and magnetic G M (q 2 ) form factors of the nucleus. G M (0) = 1 + κ, with the anomalous magnetic moment of the nucleus defined as [54] 1 + κ M N = g I 2Zm p . (S14) The dominant part of the total FS correction can therefore be calculated nuclear model independently using the Zemach radius from the literature [25]. The leading Zemach correction is δ Zemach = −0.0001911 (12). For the relativistic corrections, the following formulas were derived in Refs. [51,52]: Here, γ is the Euler constant, and R 0 = 5 3 ⟨r 2 ⟩ E . Results obtained with these formulas for a homogeneously charged sphere model for low-Z nuclei are in good agreement with tabulated results in Ref. [49]. In addition to the relativistic effects, a formula for radiative corrections to the nuclear finite size δ rad based on the exponential charge and magnetization distribution model was derived in Ref. [55]. It is with R 2 D = ⟨r 2 ⟩ E 12 . For the calculation of relativistic corrections as derived in Ref. [51], we relied on different models of the charge and magnetization distributions. In Ref. [56], parametrizations of both charge and magnetization distributions are given which are based on elastic scattering data. In Ref. [57], a charge distribution is parametrized as a sum of Bessel functions. To obtain an estimate for the nuclear model uncertainty, we also calculated the relativistic corrections using different charge and magnetization model distributions suggested for light elements in Ref. [51], like the exponential and Gaussian distributions. While using these distributions, we adjusted the model parameters such that the rms charge ⟨r 2 ⟩ E = 1.973 (14) fm and magnetic radii ⟨r 2 ⟩ M = 1.976(47) fm match the values given in Ref. [25]. In an attempt to give a conservative estimate for the upper boundary of the nuclear model uncertainty, we also used several nuclear models which are more unrealistic for low Z, that is, the homogeneous sphere and the two-parameter Fermi distribution with a skin thickness parameter of 2.3 fm [37]. For all models, relativistic corrections were similar, ranging from −2.82 × 10 −6 to −2.87 × 10 −6 . We would like to point out that the nuclear-model-dependent part of the relativistic corrections [equation (16)] amounts to less than 10%, namely, −0.22 × 10 −6 to −0.27 × 10 −6 . This total model-dependent part of the relativistic corrections is also significantly smaller than the uncertainty due to the Zemach radius. We estimated our error due to inconsistent usage of nuclear model parameters from different references as the maximum difference between any two results obtained with the above mentioned models. This error of 5 × 10 −8 is much smaller than the uncertainty due to the Zemach radius.
The uncertainty due to the rms and magnetic radii given in Ref. [25] are 5 × 10 −9 and 2 × 10 −9 , respectively. Furthermore, we estimated the uncertainty due to uncalculated higher orders in Zα to be 1.5 × 10 −8 . In total, we find the following finite-size correction to the hyperfine splitting δ FS = −0.0001939 (12) . (S18) Here, the bracket shows the uncertainty due to the Zemach radius, whereas the uncertainty originating from the difference between nuclear models is much smaller. As can be seen from equations (S12) and (S16), nuclear modeldependent quantities, i.e. the above mentioned moments of charge and magnetization distributions, are suppressed by a factor (Zα) 2 compared to the leading Zemach term, explaining the small nuclear model uncertainty. The total finite nuclear size correction corresponds to a frequency shift of 1.678 MHz, the uncertainty due to the Zemach radius being 11 kHz.
The QED correction to the hyperfine splitting is [58] δ QED = a e + δ QED,bind . (S19) Here, a e is the free-electron magnetic anomaly. Binding corrections δ QED,bind at the one-loop level were calculated to all orders in Zα using formulas and tabulated values obtained for point nuclei in Refs. [37,59]. We assumed an uncertainty of 1 in the last digit given where no uncertainty was specified in these references. At the two-loop level, binding corrections were taken into account up to O (Zα) 2 [52,58,60]. The uncertainty of the numerical coefficients in these references corresponds to an uncertainty of δ QED of 4 × 10 −9 . We estimate the uncertainty of δ QED due to uncalculated higher-order QED corrections to be 6 × 10 −9 . Muonic δ µVP and hadronic δ hadVP vacuum polarization corrections were estimated using formulas from Refs. [55,61,62]. The electroweak interaction contribution δ ew was calculated using formulas from Ref. [63].
For the calculation of the recoil correction to order O (Zα) [52,55], we use the approach developed for the hydrogen atom in Ref. [64][65][66][67]. This formalism depends on the electric and magnetic form factors of the nucleus.We calculated the recoil correction using model form factors from Ref. [56] which are based on elastic electron scattering experiments. We also used model formulas for form factors previously used in hyperfine splitting calculations for the hydrogen atom [52,68]. In this case, we adjusted the model formulas such that the rms ⟨r⟩ 2 and magnetic radii ⟨r⟩ 2 correspond to the values given in Ref. [25].  in Table 3. We assume that by using different model form factors, including the unrealistic Fourier transform of the homogeneous sphere model, we overestimate the nuclear model uncertainty rather than underestimating it.
(Especially, since the model uncertainty in Table 3 is larger than the uncertainty due to fit parameters from Ref. [56].) We calculated the uncertainties of the recoil parameter δ recoil due to the uncertainty of the magnetic radius to be 6 × 10 −8 (second uncertainty in Table S3). The uncertainty due to the charge radius is too small to be shown in Table S3. The higher-order recoil correction to order O (Zα) 2 [52] was calculated using formulas developed for the hydrogen atom in the point-nucleus approximation in Ref. [65] by replacing α → Zα. To account for uncalculated higher orders, we assign a 100% uncertainty to this contribution (third uncertainty in Table S3). We would like to point out that the nuclear model uncertainty is smaller than the uncertainty due to uncalculated higher orders.
Calculating the radiative recoil parameter using equation (2) in Ref. [62], we obtain 1.1 × 10 −7 , which corresponds to 1 kHz. Since this formula is developed for a structureless nucleus and therefore not applicable for 3 He + , we use this result as an estimate of the uncertainty due to radiative recoil (fourth uncertainty in Table S3).
Another nuclear structure contribution originates from nuclear polarization. The evaluation has been performed in the framework of a semi-analytical model, introduced in Ref. [69] and later used in Refs. [70][71][72][73]. In our recent work [74] we applied this formalism to calculations of NP corrections to the energies and g factor of H-like ions. In the case of HFS, the NP correction can be expressed analogously by the formula: .
Here, the summation goes over all nuclear excitations with energies ω L and reduced transition probabilities B(EL).
The Hamiltonian H µ = e 4π µ · r×α |r| 3 describes the interaction of the electron at the coordinate vector r with the nuclear magnetic moment µ, with the α being the usual 3-vector of Dirac matrices. n and k denote here the electronic states in the spectrum of the bound Dirac equation, the Y LM are spherical harmonics, and the radial functions F L are given in the sharp surface approximation [69] as with R 0 being the nuclear radius. The first term in the brackets corresponds to the irreducible (for k ̸ = a) and reducible (k = a) contributions, and the second one to the vertex contribution. The factor of 2 in the first term originates from the existence of two equivalent diagrams. We would like to stress that the used formalism rigorously accounts for the electronic structure and spectrum, which is important in the case of the hyperfine interaction.
The nuclear parameters ω L and B(EL) for nuclear states have been calculated in the framework of the random phase approximation (RPA) with Skyrme-type nuclear interactions [75,76]. Recent calculations [77] of the NP correction to the energy levels for different models showed the model dependence on the ≈ 15% level for heavy muonic atoms. The results are sensitive to the nuclear radius and to the use of a smoothing function F BW [78], which allows to remove the singularity of the H µ potential at the origin. The results for the relative NP corrections to the HFS are presented in Table S4. A comparison with other literature values for nuclear polarization corrections to electronic [71] and muonic energy levels -with both light [79] and heavy [80] nuclei -and to g factors and hyperfine splittings in highly charged ions [73] shows that our model yields an upper estimate of the effect in the HFS of 3 He + , and it turned out to be negligible compared to other theoretical contributions.
Analyzing available NP calculation results [79], one can see that the dominant contribution is originated from E1 transitions, fully taken into account in our approach as well. Additionally, we used full relativistic approach for the electronic wave functions. Despite the fact that RPA has been developed for heavy nuclei and should be used with caution for low-Z nuclei, we see a good agreement between our results and those available in the literature, including light nuclei [79]. This, and also the fact that the estimated NP contribution is two orders of magnitude smaller than the total theoretical uncertainty, justify the current use of our model as an approximation. Nonetheless, we conservatively assign the same large 11-kHz uncertainty to NP as to the dominant nuclear effect, namely, the finite-size effect.
With all above mentioned corrections, we obtain a total theoretical value for the hyperfine splitting of

Evaluation of the shielding constant
The combined effect of the nuclear magnetic field and the external magnetic field is expressed by the nuclear shielding contribution. The energy shift due to such diagrams was calculated in Ref. [84]. In the hydrogenlike ion, the shielding contribution is parameterized by replacing the bare nuclear g I factor by the observable, shielded nuclear g I factor, g I → g I (1 − σ), with the shielding constant σ. The factor (1 − σ) describes the modification of the nuclear magnetic moment through the presence of the bound electron.
Different contributions to the shielding constant σ were calculated following the theoretical methods in Refs. [84,85]. The leading relativistic shielding term can be given as a closed analytical expression [81]. The second largest correction comes from the nuclear recoil effect, while the dominant QED corrections are the one-loop self-energy diagrams. We found that tabulated results for the Bohr-Weisskopf correction from Ref. [84] for the ions 17 O 7+ , 43 Ca 19+ and 73 Ge 31+ can be parametrized as α(Zα) 4 (m e r rms )F BW (Zα) and found that the function F BW (Zα) only weakly depends on Z. For relatively low Z between Z = 8 and Z = 32, F BW (Zα) varies between -2.0 and -2.4 (using rms radii from [86]). Using values for F BW (Zα) in this range would correspond to a Bohr-Weisskopf correction for 3 He which is significantly smaller than our estimated uncertainty of the QED correction. We therefore assume the Bohr-Weisskopf correction to be small and give the value obtained in this way in Table S5, with a 100% uncertainty in order to give a conservative upper estimate of the possible error. In a similar way, finite-size corrections, as obtained from tabulated values in Ref. [84], can be parametrized as σ FS = α(Zα) 4 (m e r rms )F FS (Zα). Similarly to the Bohr-Weisskopf correction, we found that this function F FS (Zα) varies little between Z = 8 and Z = 32, between -2.7 and -3.9. We therefore gave an estimate for the FS correction using above formula for Z = 2 with typical values observed for F FS (Zα) which again turned out to be smaller than the uncertainty of the QED effect, and, in particular, much  smaller than the uncertainty needed to extract the nuclear magnetic moment from the experimental data. Just as in the case of the Bohr-Weisskopf correction, we indicate the estimate for the FS correction obtained in this way in Table S5, with a conservative 100% uncertainty. Attempts to model the numerical data for FS and Bohr-Weisskopf corrections using modified formulas with log(Zα) terms, α(Zα) 4 (m e r rms )F FS/BW (Zα) log(Zα), lead to functions F FS/BW (Zα) which vary much stronger than in the above case. We therefore assume that the leading order terms of these corrections do not contain logarithmic terms. Results are listed in Table S5.
A very recent calculation of the shielding constant [87] is in minor disagreement with our value, primarily due to an updated QED theory of order α 2 (Zα) 3 .

Determination of the Zemach radius
Instead of comparing experimental and theoretical values for HFS, one can alternatively determine nuclear parameters such as the Zemach radius, assuming the correctness of theory and employing the current experimental value for the nuclear g I factor. This was performed previously for the Zemach radius of the proton in Ref. [52,68]. The extraction is based on the HFS formula: Here, E F is the Fermi hyperfine splitting. This equation can be solved for the nuclear structure-dependent parameters as follows: The Zemach radius can be adjusted such that above equation is fulfilled. Given a specific set of nuclear model form factors to be used in the leading-order recoil calculation, and assuming the rms charge radius ⟨r 2 ⟩ E to be correct, one can determine both r Z and ⟨r 2 ⟩ M iteratively, starting from the current literature values. To this end, one employs the definition (S13) of the Zemach radius to determine ⟨r 2 ⟩ M from r Z . For the set of model form factors previously used in the proton Zemach radius determination [68], we find a Zemach radius of 2.608(24) fm, with the magnetic radius being 2.15(4) fm. Our extracted Zemach radius disagrees with the value of 2.528 (16), determined in Ref. [25] from electron scattering data. The uncertainty accounts for the difference of extracted Zemach radii for all different model form factors used in the leading order recoil calculation, as mentioned above, as well as the uncertainties of all remaining HFS contributions from Table S3.